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ABSTRACT 


"■  ) 

Often,  in  solving  an  elliptic  equation  with  Neumann  boundary  conditions,  a  compatibility 
condition  has  to  be  imposed  for  well-posedness.  This  condition  involves  integrals  of  the 
forcing  function. 

When  pseudospectral  Chebyshev  methods  are  used  to  discretize  the  partial  differential 
equation,  these  integrals  have  to  be  approximated  by  an  appropriate  quadrature  formula. 
The  Gauss-Chebyshev  (or  any  variant  of  it,  like  the  Gauss-Lobatto)  formula  can  not  be  used 
here  since  the  integrals  under  consideration  do  not  include  the  weight  function.  A  natural 
candidate  to  be  used  in  approximating  the  integrals  is  the  Clenshaw- Curtis  formula,  however 
we  show  in  this  paper  that  this  is  the  wrong  choice  and  it  may  lead  to  divergence  if  time 
dependent  methods  are  used  to  march  the  solution  to  steady  state. 

We  develop,  in  this  paper,  the  correct  quadrature  formula  for  these  problems.  This 
formula  takes  into  account  the  degree  of  the  polynomials  involved.  We  show  that  this 


formula  leads  to  a  well  conditioned  Chebyshev  approximation  to  the  differential  equations 


and  that  the  compatibility  condition  is  automatically  satisfied. 
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INTRODUCTION 

We  deal  here  with  a  problem  encountered  in  the  solution  via  Chebyshev  spectral  colloca¬ 
tion  discretization  of  the  Poisson  equation  with  homogeneous  Neumann  boundary  conditions. 
The  problem  arose  in  the  context  of  solution  of  the  pressure  Poisson  equation  in  a  time-split 
algorithm  for  the  incompressible  Navier-Stokes  equations.  For  this  problem  to  be  well-posed, 
the  source  term  must  satisfy  a  compatibility  condition;  the  numerical  analog  of  this  condi¬ 
tion  using  straightforward  Clenshaw- Curtis  quadrature  formulae  was  found  to  be  numerically 
ill-conditioned  and  produced  large  distortions  in  the  spectral  solution. 

In  Section  I,  we  outline  the  time-split  algorithm  and  describe  the  difficulty  encountered; 
we  proceed  in  Sections  II  and  III  to  analyze  the  difficulty,  first  in  terms  of  an  equivalent 
parabolic  equation  which  would  be  utilized,  for  instance,  in  the  iterative  solution  of  the 
Poisson  equation  of  interest,  then  in  terms  of  the  steady  equation  itself.  The  proper  quadra¬ 
ture  formula  is  also  developed,  which  alleviates  the  numerical  difficulty.  In  Section  IV  we 
show  numerical  examples  to  demonstrate  the  numerical  inconsistency  and  its  resolution. 

I.  TIME-SPLIT  ALGORITHM 

A  splitting  method  is  employed  in  many  simulations  to  advance  the  solution  of  the  incom¬ 
pressible  Navier-Stokes  equations  from  time  t"  to  Writing  the  Navier-Stokes  equations 
in  vector  notation,  with  u  representing  the  velocity  {u,v,w)  we  have: 

Ut  +  u  •  Vu  =  —  VP -t- i/V^u  (l.lo) 


and 

Vu  =  0  (1.16) 

in  the  domain  fl,  and  : 

u  =  0  on  the  boundary  P 

The  split  scheme  first  advances  u"  to  an  intermediate  solution  u*  by  solving: 

u; -f  u*  •  Vu*  =  i/V^u*  (1.2) 

u*  =  g*  on  r 

The  intermediate  boundary  condition  u*  =  g*  is  discussed  in  [1].  Finally,  the  solution  is 
advanced  from  u*  to  via: 

u;*+i  =  _vP"+‘ 

Vu"+^==0  (1.3) 

n  -  =  0  on  r 
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where  n  is  the  unit  normal  to  the  boundary  F.  Note  that  the  final,  “pressure  correction” 
step  by  itself  is  a  set  of  in  viscid  equations;  and  is  well-posed  when  boundary  conditions  on 
the  normal  component  of  the  velocity  only  are  enforced.  At  the  end  of  the  full  step  there 
exists  a  non  zero  tangential  component  of  velocity  on  the  boundary.  The  magnitude  of  this 
slip  velocity  „an  be  reduced  to  O(Ai^)  by  a  proper  choice  of  the  intermediate  boundary 
condition  on  u*  [1]. 

Using  backward  Euler  time  discretization  for  Eq.  (1.3)  yields: 

=  u*  -  (1.4) 

The  pressure  step  is  actually  carried  out  in  two  parts.  First,  the  divergence  of  Eq.  (1.4) 


yields; 


^2pn+l  ^  y  . 

At 


where  V  •  =  0  is  enforced.  Then  the  velocities  are  updated  using  Eq.  (1.4). 

Note  that  this  formulation  requires  a  boundary  condition  for  the  pressure.  This  poses  a 
problem,  since  there  is  no  natural  boundary  condition  for  pressure.  The  use  of  a  condition 
derived  from  enforcing  the  normal  momentum  equation  at  the  boundary  was  attempted,  but 
was  apparently  inconsistent  as  it  resulted  in  explosive  instability.  Fortunately,  analysis  of 
the  split  scheme  itself  yields  a  self-consistent  pressure  condition: 

n-VP"+Mr=0  (1.6) 

since  both  n-  u*  and  n-  are  zero  on  the  boundary.  The  error  involved  in  this  specification 
is,  we  believe,  related  to  the  overall  splitting  error  of  the  scheme.  It  is  also  known  that  the 
error  due  to  imposition  of  an  inconsistent  Neumann  pressure  boundary  condition  is  isolated 
to  a  thin  “boundary  layer”  [4]  . 

It  is  easily  shown  using  the  divergence  theorem  that 


^(V  u*)dn  -0 


is  required  for  the  pressure  Poisson  equation  to  be  well-posed  with  this  boundary  condition. 
For  application  of  this  algorithm  in  closed-system  problems  with  homogeneous  Dirichlet 
velocity  boundary  conditions,  numerical  tests  indicate  Eq.  (1.7)  is  satisfied  in  general.  How¬ 
ever,  in  applications  of  a  recently- developed  non-reflecting  outflow  boundary  treatment  [5] 
to  simulations  of  flow-through  systems,  it  was  found  that  Eq.  (1.7)  did  not  hold  discretely 
for  the  intermediate  velocity  field.  As  will  be  shown  in  the  numerical  examples  of  Section 
IV,  this  leads  to  large  distortions  in  the  computed  pressure  field.  It  was  decided  that  since 
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there  is  no  solution  to  the  problem: 


V^P  =  cast  7^  0  in  n 
V  P  •  n  =  0  on  r 


the  pressure  Poisson  equation  would  be  modified  to  read: 


V^P  = 


(1.8) 


The  integral  in  Eq.  (1.8)  was  implemented  using  the  straightforward  Clenshaw-Curtis 
quadrature  formulae,  appropriate  for  Chebyshev  collocation.  Although  the  pressure  field 
distortions  were  much  reduced,  their  magnitude  was  still  unacceptable.  The  analysis  pre¬ 
sented  in  the  following  show  the  inconsistency  in  the  use  of  these  formulae,  and  the  proper 
quadrature  formula  to  recover  a  smooth  solution  to  Eq.  (1.8). 
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II.  TIME  DEPENDENT  PROBLEMS 


Consider  the  parabolic  equation 


Ut  g 


(2.1a) 


u(x, 0)  =  0 

with  the  Neumann  boundary  conditions 

u*(l,t)  =  =  0. 


(2.16) 


(2.2) 


We  assume  also  that 


so  that 


J ^ g(x)dx  = 


u(i,  f)<fx  =  tXi(l)  -  Ux(-l)  +  J  ^dx  =  0 
and  therefore  in  view  of  (2.1b) 

J  u(x,t)dx  =  0  for  all  t. 


(2.3) 

(2.4) 


(2.5) 


In  the  pseudospectral  Chebyshev  method  we  seek  an  i  polynomial  u/^(x,t),  that  satisfies 
the  boundary  conditions  (2.2),  such  that  (2.1a)  is  satisfied  at  the  points  Xj  =  cos  namely 
we  seek  ut/{x,t)  such  that 


duff 

~w 


■5i5-  +  s 


at 


X  —  Xj  j  —  .  ,N  —  \ 


UAr(x,0)  =  0 


(2.6a) 

(2.66) 


and 

^(_1.0=^(M)  =  0.  (2.7) 

We  refer  the  reader  to  [3]  for  a  discussion  of  the  stability  and  convergence  of  (2.6)  and  (2.7) 
to  (2.1)  and  (2.2). 

Here  we  are  interested  in  the  question  whether  the  numerical  approximation  satisfies 
the  conservation  property  (2.5).  Since  the  numerical  solution  Uff{x,t)  defined  in  (2.6)  is  a 
polynomial  of  degree  in  x  it  is  natural  to  consider  the  Clenshaw  Curtis  quadrature  formula 
PI-  It  uses  the  collocation  points  of  (2.6a)  and  it  is  exact  for  polynomials  of  degree  N. 
Lemma  (1.1).  Let  N  be  an  even  number  and  let  qj  be  defined  by 


a 


4  ’  1 


k=0 


2iTlk  1  „  , 


(2.8) 
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Oio  —  oti^  —  — -  where  Cq  =  Cff  =  2  Cfc  =  l  Q  ^  k  ^  N 

—  1 


Then  for  any  polynomial  h{x)  of  degree  at  most  N 

/  h{x)dx  =  h(xi)ai 


where  xj  =  cos  ^ . 

A  close  inspection  of  (2.6a)  reveals  that,  in  general,  ujv(x,t)  does  not  satisfy  the  con¬ 
servation  condition  (2.5).  Indeed,  since  g(x)  in  (2.1)  is  not  a  polynomial  of  degree  N,  the 
quadrature  formula  (2.9)  is  not  exact  for  g,  whereas  it  is  exact  for  the  other  terms  in  the 
equation. 

To  remedy  the  situation  it  is  customary  to  modify  (2.6a)  by  seeking  Uf;(x,t)  that  satisfies 


du^  d^Uf4  1 


^  +  at  x  =  x,  l<j<N 


(2.6c) 


where  q;  are  defined  in  (2.8).  Of  course  U{i{x,t)  still  satisfies  (2.1b)  and  (2.2).  Thus  the 
right  hand  side  of  (2.6c)  has  zero  mean  if  the  Clenshaw- Curtis  formula  is  being  used.  Unfor¬ 
tunately,  even  with  the  above  modification  the  solution  of  (2.6c)  does  not  satisfy  the  discrete 
analog  of  (2.5);  in  fact  we  can  state; 

Lemma  (1.2):  Let  u^{x,t)  be  the  solution  of  (2.6c,b)  with  the  Neumann  boundary  condi¬ 
tion  (2.7).  Let  gN{i)  be  the  coefficient  of  the  TV^’th  Chebyshev  polynomial  in  the  expansion 
of  g{x,  i),  namely 


Then 


where 


-  1  yi  ff(x)Tiv(x)j_ 

9n  -  /  /r - 5-  ax 

TT  J-i  vl  -  X^ 

(2.10) 

T/v(x)  =  cos  N(cos~^  x). 

(2.11) 

t)aj  =  ^2  _  1  {“^{0  1  9N{t)dT'^ 

(2.12) 

%(()  -  ‘  /‘ 

TT  J-l 

(2.13) 

Remark:  Equation  (2.12)  points  out  the  difficulty  in  the  method  (2.6c).  Suppose  for  exam¬ 
ple  that  ^  is  a  time  independent  function,  then  the  solution  of  (2.1a,b)  converges  to  a  steady 
state  solution  and  therefore  its  highest  coefficient  un(0  is  a  bounded  function  of  t.  However 
the  second  term  in  (2.12)  is  not;  in  fact 


I  gN{r)dT  = 
Jo 
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and  this  diverges  linearily  in  t. 

The  above  does  not  contradict  the  usual  stability  argument  since  for  fixed  t,  UN{x,t) 
converges  to  u(x,t)  as  N  tends  to  infinity  ,  but  the  usefulness  of  (2.6c)  as  a  method  for 
marching  to  steady  state  is  doubtful.  This  will  be  demonstrated  further  in  the  numerical 
cxpcnrp-cn^s. 

Proof;  We  denote  by  Ptjg{x,t)  the  (unique)  interpolation  polynomial  of  degree  iV  in  x  that 
interpolates  g{x,t)  at  the  points  i/  =  cos  0  <  /  <  iV.  With  this  notation  equations 
(2.6c)-(2.7  )  can  be  rewritten  as 

~W  ^  ^  +  {'4“'+  B)T'„(x)  (2.14o) 

»;v(x,0)  =  0  =  ^(1,0  =  (2.144) 

By  comparing  the  coefficients  of  the  highest  Chebyshev  polynomial  in  both  sides  of  (2.14a) 
one  gets 

A  = 

Multiplying  now  (2.14a)  by  and  summing  we  get 


duff 

dt 


-  9N 


N' 


(2.16) 


d 

—  yuN{Xj,t)a, 

j~0 


N 

E 

j-ij 


UN 


N 


dx^ 


{x„t]aj  +  '^g(xj)aj 


J^g(x()ai  +  '^(Axj  +  B)T'yxj)aj 


1=0  j=0 

(2.16) 


We  use  now'  the  exactness  of  the  Clenshaw  Curtis  formula  and  the  fact  that  T'yx,)  =  0  1  < 

j  <  A'’  —  1  to  get 


d  duN,  . 


du 


N 


dx 


(  — l,t)  +  2AN^ao 


2iV2 
-  1 


diiff 


-  9s{t) 


1 

N 


(2.17) 


Integration  of  (2.17)  yields  (2.12)  and  the  proof  is  completed. 

It  is  obvious  that  the  problem  lies  in  the  use  of  the  Clenshaw  Curtis  quadrature  formula. 
Specifically,  we  need  to  use  an  accurate  quadrature  formula  that  does  not  use  the  boundary 
points.  Such  a  formula  is  constructed  in  Lemma(1.3). 

Lemma  (1.3):  Let  /(x)  be  a  polynomial  of  degree  JV  —  2  at  most.  Let  Xj  =  cos  ^  then 


/  I{x)dx  =:  g(x;)/5j 


(2.18a) 


where 
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[2.184) 


and  aj  are  defined  in  (2.8). 

Proof:  Since  /(x)  is  a  polynomial  of  degree  iV  —  2  ,  it  is  uniquely  determined  by  its  values 
at  the  interior  points  Xj,j  =  1,. . .  ,N  —  1.  In  fact 


Therefore 


where 


j=l  ^  ^3 


I  f{x)dx  =  Y1  fi^3)0j 

,=-.1 

fl  t'  tx)  1 

P3  =  /  •  7f;^Jx. 

J~\x  —  Xj 


We  use  now  the  Clenshaw  Curtis  formula  to  evaluate  the  integral  (2.20).  Thus 

~  ^3  Tn(^])  ■’ 

Mi 


(2.19) 


(2.20) 


T'M  TUz}h^ 

(i-x^mx,)  °-"(-i-x,)r;^(x,) 


(2.21) 


=  (-ly 


and  the  quadrature  formula  (2.18)  is  thus  established. 

Using  Lemma  (1.3)  we  suggest  the  following  algorithm  for  the  pseudospectral  Chebyshev 
discretization  of  (2.1)-(2.2).  In  the  new  algorithm  we  seek  ui^{x,t)  such  that 


duN  1  f  ^  -1  M  ^ 

~dr  ^  ~d^  +  5  “  b  L  9{^3P3  at  X  =  X,  j  =  1, . . . ,  iV  -  1 


(2.22a) 


(2.226) 


ujv(x,0)  =  0.  (2.2c) 

In  the  new  method  the  compatibility  condition  (2.5)  is  satisfied.  In  fact  we  state 
Theorem  (1.1):  Let  ui^{x,t)  be  the  solution  of  (2.22).  Let  be  defined  in  (2.21)  then 
^N(xj,t)/3j  =  0  for  all  t. 
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Proof;  We  multiply  (2.22a)  by  and  sum  to  get 


j  N-l  N-1  ffl.  N-\  N-\ 

—  Y,  UN{xj,t)/3j  -  Y  H  -  Y1  9{^3)Pr  (2-23) 

j=i  j=i  j-i  j=i 


We  use  now  Lemma  (1.3)  and  the  fact  that  is  a  polynomial  of  degree  N  —  2,  to  get 


/•!  d^UN  Oun,  .  Oun  ,  . 


(2.24) 


and  therefore 


and  since  UAr(x,0)  =  0 


which  proves  the  theorem. 


d 

^  ^  '^Ni^jyi)0j  —  0 
N-l 

Y  =  0 

i=i 


(2.25) 
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III.  STEADY  STATE  PROBLEMS 

In  this  section  we  discuss  the  equation 


—  Q 


(3.1a) 


7i^(±l)  =  0  (3.16) 

Equation  (3.1)  can  be  viewed  as  the  steady  state  version  of  (2.1a),  however  this  time  we 
need  one  more  condition  to  get  a  unique  solution.  We  impose  the  condition 


j  u  dx  =  0. 


In  order  for  (3.1a)  and  (3.1b)  to  be  compatible  g[x)  has  to  satisfy  the  condition 


J  g(x)dx  —  0. 


We  have  demonstrated  in  Section  II  that  the  use  of  the  Clenshaw- Curtis  quadrature 
formula  (2.8)  causes  problems  ii.  solving  (3.1)  by  trying  to  evolve  the  solution  of  (2.1)  to 
steady  state.  Here  we  would  like  to  discuss  the  direct  solution  of  (3.1)  via  the  influence 
matrix  technique.  We  will  demonstrate  that  the  quadrature  formula  (2.18)  developed  in 
Section  II  should  be  used,  rather  than  the  Clenshaw- Curt  is  formula  (2.8). 

The  solution  of  (3.1)-(3.2)  via  the  influence  matrix  technique  involves  seeking  the  ap¬ 
proximation  of  (3.1)  as  a  sum  of  the  solution  of  three  problems: 

Problem  1:  Seek  a  polynomial  W[^{x)  such  that 

— TT"  =  9  -  9o  at  X  =  ij  =  cos  ■—  1  <  J  <  iv 
dx^  iV 

u;iv(l)  =  0  =  ii;Ar(-l)  (3.4) 

and  go  is  an  approximation  to  / g[x)dx.  We  will  consider  either 


9o  = 

j=o 


Uj  given  in  (2.8) 


(3.5a) 


N-l 

9o=  9{xjWj- 

f3j  given  in  (2.18b) 

Problem  2:  Seek  a  polynomial  v^{x)  such  that 


UtLf  m  r  ^ 

=  0  at  X  =  Xj  J  =  1, . . . ,  N  —  1 
dx^ 


(3.56) 


(3.6a) 
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=  0  v^(l)  =--  1 

Problem  3;  Seek  a  polynomial  v^^(x)  such  that 


Y-  =  0  at  X  =  Xj  j  =  1, . . .  ^  JV  -  1 


(3.6b) 


(3.7a) 


=  1  v^^(l)  =  0 

Given  Wfj we  look  for  a  solution  of  the  form 


(3.76) 


un(x)  =  wn{x)  +  div'{x)  +  d2V^\x). 


Clearly  uw  satisfies 


d^U[^(x) 


=  g(x)  —  gc  at  x  =  Xj  I  <  j  <  N 


for  any  dijd^.  The  constants  di,d2  are  determined  by  the  boundary  conditions  (3.1b). 
In  our  case  (3.1b)  (3.2)  imply 


n  <iuN(l)  dwt^{l)  _  dv^(l) 

~  j  ^N{^)dx  +  di  j  v\x)dx  +  dj  j  v^\x)d3 


(6)(3.9) 


Ideally  (3.9a)  and  (3.9b)  should  yield  the  same  equation  for  di,  dj.  The  integrals  in  (3.9c) 
have,  or  course,  to  be  discretized. 

In  this  simple  model  problem  we  are  able  to  evaluate  explicitly  Wf^(x),v^ (x),  and  v^^(x). 
In  fact  we  can  state 

Lemma  (2.1):  The  solutions  of  problems  2  and  3  are  given  by 


/  1  +  X  jj 

V  =  -  V 

2 


1  —  I 


(3.10) 


Proof:  i;^(x),  u^^(x)  in  (3.10)  are  the  only  polynomials  that  satisfy  (3.6),  (3.7)  respectively. 
In  order  to  get  an  expression  for  wi^(x)  we  introduce  the  interpolation  polynomial  Fn-25 
that  interpolates  g(x)  at  the  points  Xj  1  <  j  <  iV  -  1.  By  (2.19) 


,  3  'ic''  ,  1 

9(^)  =  E  ■  rj^Y 


(3.11) 
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It  is  clear  form  (3.-1)  that  satisfies 


dx^ 


-29  go 


(3.12a) 


ii'w(  1 ) -- 0  -  I )-  (3.126) 

Note  that  equation  (3.12a)  holds  for  any  -I  <  x  <  1,  not  just  for  the  grid  points  Xj\  it  is 
therefore  possible  to  integrate  (3.r2al  to  get 

^  J ~ 

and 

(ixu  hi 

Wn{x)  -  -^(  -1)(1  i  i  y  ™  7/)[/^Ar_o^(7/)  --  ^0^7/.  (^-1^0 

Substituting  now  the  condition  u;yv(l)  --  0  we  get 

^^(-1)  ~\J_y  ~v)y'N  -2g{v)  -  go]dri.  (3.15) 


We  are  ready  now  for  the  main  result  of  this  section. 

lA'inma  (2.2):  Let  go  be  defined  by  the  use  of  the  Clcnshaw-Curtis  formula  as  in  (3.5a), 
then  e(|uation  (3.9a)  is  incompatible  with  (3.9b). 

Proof:  Using  the  expressions  (3.10)  for  and  in  (3.9a, b)  we  get 


dwf^ 

dx 


(-0 


dj_ 

2 


d2 

T 


However  from  (3.13) 


dtUN 

dx 


(1) 


di  c/2 

~2  Y' 


dwt4 

dx 


(0 


^(-1)3  /V^_2.9(0-5oK- 


(a) 

(3.16) 

{h) 


(3.17) 


Since  Pn  29  's  a  polynomial  of  degree  N  -  2,  the  quadrature  formula  (2.l8)  is  exact  and 
therefore  ^ 

j  y^N  2.g(0  -  go\di  Y,g{-^j)Pj  -  (3-i8) 


j  1 


N 


g{xo)oo  I  g{xf])(^N 


j  o 


2^ 


j  ‘> 


11 


where  cq  =  2  =  c^,  Cj  —  1 
so  that 

d'^'N  ,  K  .  dwN  , 

which  proves  the  Lemma. 

However  if  one  uses  the  new'  quadrature  formula  as  in  (3.5b)  one  gets 
Lemma  (2.3):  Let  go  be  defined  in  (3.5a),  then  (3.9)(a)  and  (3.9)(b)  are  compatible. 
Proof;  We  start  as  in  the  last  Lemma  except  now 

/  {Pn-29  -  go)d(  =  9{xj)(3j  -  =  0 

which  yields  the  result. 

We  have  thus  established  that  the  use  of  the  Clenshaw- Curtis  formula  may  lead  to  an 
approximation  that  does  not  satisfy  (3.1b).  The  new  quadrature  formula  alleviates  this 
problem. 
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IV.  NUMERICAL  RESULTS 

In  Section  II  we  showed  that  the  use  of  the  Clenshaw-Curtis  formula  (2. 8), (2. 9)  may  cause 
problems  when  one  attempts  to  solve  the  steady  state  problem  corresponding  to  (2.1).  In  this 
Section  we  report  on  the  numerical  solution  of  equation  (2.1),  (2.2)  with  two  different  source 
functions,  g(x)  —  cosStti  and  g{x)  =  cos  5x1.  To  advance  in  time  we  use  the  fourth  order 
Runge-Kutta  scheme;  a  grid  of  19  points  was  used.  In  table  I  we  summarize  the  results  for 
the  first  case:  g[x)  =  cos  3xx.  The  first  column  gives  the  L2  deviation  from  the  steady  state 
solution  when  (2.6a,b)-(2.7)  is  used,  i.e.  no  modification  for  the  right  hand  side.  A  linear 
growth  is  observed  m  the  deviation  from  the  steady  state.  The  solution  did  not  converge 
even  after  200,000  time  steps. 

In  the  second  column  we  give  the  results  for  (2.6c),  in  which  the  approximation  to  the 
mean  of  g  is  obtained  with  the  use  of  the  Clenshaw-Curtiss  formula.  It  is  interesting  to  note 
that  the  results  are  essentially  the  same;  the  subtraction  did  not  improve  the  situation.  In 
the  third  column  we  present  the  results  for  (2.22a);  convergence  was  obtained  after  10,000 
time  steps. 

In  table  II  we  present  the  same  results  for  g(x)  =  cosSxx.  Basically,  the  results  are  the 
same  as  in  table  I;  however,  the  divergence  is  much  more  rapid  than  in  the  first  case.  The 
new  method  converged  after  15,000  time  steps. 
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Table  I 

L2  Errors  for  g(x)  =  cosStti  N  —  19 


No.  Of 

No 

Clenshaw- Curtiss 

GS 

Time  Steps 

1  Subtraction  | 

Subtraction 

Subtraction 

1.778 

(-3) 

1.778 

(-3) 

1.778 

(-3) 

1.589 

(-6) 

8.644 

(-7) 

2.669 

(-7) 

15000 

2.324 

(-6) 

2.32 

(-6) 

2.669 

(-7)* 

1.492 

(-5) 

1.494 

(-5) 

same 

2.234 

(-5) 

2.237 

(-5) 

same 

2.3013 

(-5) 

2.9056 

(-5) 

same 

Growth 

7.4165 

(-7) 

7.4277 

(-7) 

0 

converges  after  25000  steps 


Table  II 

Z/2  Errors  for  g{x)  =  cosSttx  =  19 


No.  Of 

1  No 

Clenshaw-Curtiss 

1  GS 

Time  Steps 

Subtraction 

Subtraction 

1  Subtraction 

6.320 

(-4) 

6.320 

(-4) 

6.320 

(-4) 

1.4509 

(-3) 

1.456 

(-3) 

2.4405 

(-4) 

15000 

2.1173 

(-3) 

2.126 

(-3) 

2.4554 

(-4) 

100000 

1.3511 

(-2) 

1.3606 

(-2) 

same 

150000 

2 

(-2) 

2.036 

(-2) 

same 

200000 

2.581 

(-2) 

2.644 

(-2) 

same 

Growth 

6.46 

(A) 

6.75 

r-4) 

Finally, 

we  show  the  results  of  a  two-dimensional  calculation  for; 

'^xx  +  Uyy  =  g{x,y)  =  cos  TTX  ■  cos  Tty 

(4.1a) 

with 

Uj:  =  0  on  (±l,y) 

(4.16) 

Uy  =  0  on  (x,  ±1). 

(4.1c) 

Solution  is  carried  out  by  the  influence  matrix  technique,  with  the  related  Dirichlet  Poisson 
problems  solved  via  the  direct  solution  method  of  tensor-product  diagonalization.  In  order 
to  set  the  level  of  the  solution,  that  is,  to  set  the  arbitrary  constant  in  the  solution  of  the 
pure  Neumann  problem,  one  must  set  the  value  of  the  solution  at  one  point  on  the  boundary 
by  replacing  one  row  of  the  influence  matrix.  However,  if  the  source  function  does  not  satisfy 
the  discrete  version  of  the  compatibility  relation 

J  ^  j  ^g{x,y)dxdy  =  0  (4.2) 

then  large  distortions  appear  in  the  solution  near  that  point.  On  the  other  hand,  if  (4.2)  is 
discretely  satisfled,  then  the  solution  near  the  set  point  will  be  regular. 

In  Figure  1  is  shown  isolines  of  the  solution  to  (4.1)  ,  for  which  g{x,y)  was  replaced  by 

-^J29{^hyk)akai  (4.3) 

/=o  fe=o 

that  is,  using  the  Clenshaw-Curtiss  formula  to  satisfy  compatibility.  The  location  of  the 
boundary  point  which  is  set  in  the  influence  matrix  is  obvious;  the  true  solution  is  swamped 
by  the  artifact  distortions.  In  Figure  2  is  shown  the  isolines  for  the  solution  wherein  the 
compatibility  relation  is  satisfied  using  the  quadrature  formula  developed  here.  The  solution 
is  now  smooth,  symmetric,  and  regular  near  the  set  point. 
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